Inverse design enables large-scale high-performance meta-optics reshaping virtual reality

Meta-optics has achieved major breakthroughs in the past decade; however, conventional forward design faces challenges as functionality complexity and device size scale up. Inverse design aims at optimizing meta-optics design but has been currently limited by expensive brute-force numerical solvers to small devices, which are also difficult to realize experimentally. Here, we present a general inverse-design framework for aperiodic large-scale (20k × 20k λ2) complex meta-optics in three dimensions, which alleviates computational cost for both simulation and optimization via a fast approximate solver and an adjoint method, respectively. Our framework naturally accounts for fabrication constraints via a surrogate model. In experiments, we demonstrate aberration-corrected metalenses working in the visible with high numerical aperture, poly-chromatic focusing, and large diameter up to the centimeter scale. Such large-scale meta-optics opens a new paradigm for applications, and we demonstrate its potential for future virtual-reality platforms by using a meta-eyepiece and a laser back-illuminated micro-Liquid Crystal Display.

The electric field at the target focal spot is computed using a convolution between the local sources and the appropriate Green's function �⃗ ��⃗ � = � ⃗ (�⃗, � �⃗) ⊙ ⃡ ��⃗, �⃗ � �⃗ The support of the integral is technically infinite, but a "windowing" method truncates this integral accurately to a finite region. 6 This integral is two-dimensional, which is expensive to evaluate for large area metasurfaces (diameter > 1000λ). To improve efficiency, we impose cylindrical symmetry in the parameter function, which translates into a radial local field and a radial source function. Even though the meta-atom parameters are chosen to be rotationally symmetrical, the surface is not axisymmetric at a subwavelength scale, in contrast to recent work on axisymmetric metasurfaces. 7 In our locally periodic approximation, the Green's function is also cylindrical symmetric when the focal spot is on the focal axis, so the problem is reduced to the one-dimensional convolution, which dramatically speeds up the evaluation of the intensity ( ) = � �⃗ � �� in the far field. A surrogate model-a data-driven fit that predicts the solution of Maxwell's equations instead of solving the equations numerically during optimization-is used to evaluate the local field of a given set of parameters efficiently (10 6 × speed gain, see "Surrogate model vs. RCWA" section below). Here, the surrogate model is based on Chebyshev interpolation.

Equivalence principle and Green's function for polarization conversion
is constant when is on the cylinder axis as we swipe the rotation angle of the polar coordinate and fix its radius, so the Green's function is cylindrical symmetric when the focal spot is on the focal axis of the lens. The contribution of a given radius r to the far- ( ) (� + �) where ( , ) is the position of the source in polar coordinates, which shows that the convolution boils down to a single radial term in the integrand. The reasoning is similar for a left circular polarization.

Wirtinger calculus and gradient of the intensity at the focal spot
Wirtinger calculus 8 is a common framework to differentiate complex-valued functions. Instead of taking the partial derivative of a complex number = + with respect to its real and imaginary parts, Wirtinger calculus takes the partial derivatives with respect to z and its complex conjugate ̅ = − . This formalism is especially useful for non-holomorphic functions of the complex plane like ( ) = | | 2 = ̅ , where the differential becomes = + In the case of the gradient of the intensity, the relevant complex number is the electric field E and the intensity is ( ) = ( ) � ( ), where p is the vector of all the parameters defining the metasurface.
The only thing left to do is to find the gradient ∇ ( ), starting from Equation. 1 in the main text We can transfer ∇ inside the integral sign using the dominated convergence theorem Since the Green's function is independent of the parameters of the metasurface, it can be factored out of the differentiation.
For completeness, ∇ ⃗ ( , � �⃗) is obtained by taking the analytical derivative of the surrogate model for the local field, i.e., the derivative of our Chebyshev interpolation of the zeroth order Fourier coefficient of the transmitted field with respect to the meta-atom parameter. Therefore, the surrogate model enables us to efficiently evaluate both the local field and its gradient. In contrast, methods that are not based on a surrogate model need to solve for Maxwell's equation both for the forward and the adjoint problem.

Cylindrical symmetry and Cartesian grid
To further speed up the inverse design, we use a solver with cylindrical symmetry. The parameter function we use is cylindrical symmetric; however, the evaluation on the parameter function is on a Cartesian grid. Figure S1. When defining the parameters on a Cartesian grid (blue dots) using a cylindrically symmetry function, the number of degrees of freedom (orange dots) is significantly bigger than the Cartesian grid resolution (green lines). Figure S1 shows the schematic of a metalens, which has meta-atom (blue dots) arranged in Cartesian coordinates (green grid). The spacing between the grid is the unit cell period. Each meta-atom has an "image" on radial axis (X) as represented by an orange dot. When imposing cylindrical symmetry, we consider the contribute from each meta-atom on the Cartesian grid via its "image" along the radial direction. In simulation, we have two additional sampling points within each unit cell period shown as red dots. The EM responses of the "image" meta-atoms (orange dots) are estimated through interpolation using the fine grid. In contrast to a solver, which provides the full electric and magnetic fields when given a discretized geometry, a surrogate solver returns the complex transmission when given a parameterization of the geometry. To train our surrogate model, the data was rigorously simulated on a grid with the width and the length ranging from 50 to 340 by increment of 5 nm using an RCWA method based on reticolo, 9 which results in 3481 simulations. Chebyshev regression fits a surrogate with 250 polynomials of degree up to 22, using the Chebyshev basis. The implementation uses ApproxFun.jl in Julia language, and the regression is performed by constructing the Chebyshev-Vandermonde matrix and performing a least-squares fit. Chebyshev regression smooths the data, especially where there are resonances, to reduce rapid oscillations in the surrogate model (which may result in poor local optima). The error introduced by this smoothing is negligible compared to the error of LPA or fabrication error, because the optimal design largely employs parameters far from such resonances. Evaluation of the surrogate model takes 8.5 µs on average, whereas solving with the RCWA library takes 10 s on average, which corresponds to a 6 orders-of-magnitude speedup for our surrogate-based model.

Surrogate model vs. RCWA
The experimental validation matches numerical results because the focal lengths measured correspond to the ones optimized by the inverse-design framework. This shows that LPA is accurate in the regime of our lenses' designs. Various computational methods to improve upon the accuracy of LPA exist, including overlapping-domain decomposition, 10 or spatially truncated multipole expansions 11 (which unfortunately cannot currently incorporate substrates). However, in the regime where LPA works, as in the case of this paper, the more accurate solver gives very similar results to LPA but at a major sacrifice in performance (especially when surrogates are used in LPA). LPA is quite robust to differences in neighboring atoms, and is observed to work well even for large variations in localized portions of the metasurface-it only breaks down in situations very different from those in this paper involving rapid variations over large areas, such as occur in designs for very oblique angles of incidence 12 or for very high numerical apertures. 13 Optical setup for metalens measurements Figure S3. Schematic of the optical setup used to measure the focusing performance of the metalens in three dimensions. The laser output (either from laser diodes or a super continuum laser) is collimated by a fiber coupled collimator (Thorlabs, RC12APC-P01) and incident on the metalens. The three-dimensional intensity distribution of the focal spot is then imaged through a home-built microscopic system as shown in the dashed box. It consists of a 100X objective lens (Olympus), a tube lens (Thorlabs TTL180-A), and a CMOS camera (Thorlabs, DCC1545M). The microscopic system is mounted on a translational stage. Figure S4. Schematic of the optical setup used to characterize the imaging performance of the metalens. The laser source is either a laser diode or a super continuum laser. The output light is collimated by a fiber coupled collimator (Thorlabs, RC12APC-P01) and goes through a diffuser. It is then slightly focused by an objective lens (Mitutoyo objective, 10x magnification) onto the United States Air Force resolution target. The metalens is placed at focal length away from the surface of the target. The relayed image by the metalens and the tube lens (Thorlabs TTL180-A) is captured by a CMOS camera (PCO, panda 4.2). The metalens, the tube lens, and the camera are mounted on a translational stage. Figure S5. Measured focal intensity distribution in the XZ cross section at design wavelengths when (A) incident light is in the LCP state, and output light is selected in the RCP state; (B) incident light is in the RCP state, and output light is selected in the LCP state; (C) no particular polarization state is selected for incident and output light.
We measured the focusing intensity distribution of the 2-mm-diameter, NA = 0.7, RGBachromatic metalens in three polarization configurations. We put a polarizer pair consisting of a wire grid polarizer and an achromatic quarter waveplate between the collimator and metalens to select the polarization state of incident beam, and we put another analyzer pair between the tube lens and the CMOS camera to select the polarization state of output beam. Figure S5(A) shows the measured focusing intensity distribution in the XZ cross section when we configure the incident light in the LCP state and output light in the RCP state. Figure S5(B) shows the measurement result when we configure the incident light in the RCP state and output light in the LCP state. The measurement results in two configurations are identical. It is because each metaatom has the same response when converting light from an LCP state to an RCP state or vice versa according to Jones' matrix. Since any polarization state can be decomposed into the orthogonal basis of LCP and RCP states, the metalens can work for incidence in an arbitrary polarization state and has polarization-insensitive focusing performance. Figure S5(C) (as also shown in the main text) shows the measurement results when we did not particularly select the polarization state of incident or output light (no polarizer and analyzer is used), and it shows the same focusing performance as Fig. S5(A) and (B). Figure S6. Measured focal intensity line profile (solid line (A) to (C)) at 488 nm, 532 nm, and 658 nm, respectively, in comparison with ideal Airy function profiles (dashed line). Figure S6 shows the focal intensity line profile of the RGB-achromatic metalens (NA = 0.7) at blue, green, and red wavelengths. The theoretical Airy function profiles are also shown for comparison. The calculated Strehl ratios are 0.97, 0.96, and 0.94, respectively. Here, the light sources used are three laser diodes. In comparison, we also used a super continuum laser as the light source to measure the focal spots of the metalens. Figure S7, (A) to (C), shows the focal intensity distribution in the XY cross section at design wavelengths. Figure S7, (D) to (F), shows the corresponding focal intensity line distribution. The measurement shows slightly different from the previous measurement when the laser diodes are used. There are some background light around focal spots, and the PSF profile shows focusing intensity tail decaying as radial distance in contrast with Airy function profile, which shows negligible focusing intensity within side lobes. The measured FWHMs at wavelengths of 488 nm, 532 nm, and 658 nm are 450 nm, 468 nm, and 520 nm, respectively as shown in Fig. S8. They are slightly larger than the theoretical values of Airy functions, which are 358 nm, 391 nm, and 483 nm, respectively. Figure S9 (A). Measured spectra of blue, green, and red laser diodes. (B). Measured spectra of a super continuum laser when the peak intensity of output light is tuned to center at 488 nm, 532 nm, and 658 nm (corresponding to operation wavelengths of laser diodes).

Point spread functions (PSFs) of the RGB-achromatic metalens
The measurement difference between two experiments can be explained by the laser spectra. Figure S9(A) shows the spectra of the laser diodes operating at blue, green, and red wavelength. In comparison, Figure S9(B) shows the spectra of the super continuum laser when the wavelength of the peak intensity is tuned to match with laser diodes. The laser diodes have a linewidth of 0.5 -1 nm whereas the super continuum laser has a minimum linewidth of 4 -7 nm, which is about one order of magnitude larger. The larger beam linewidth of the super continuum laser causes the broadening of the focal spots and the showing up of the focal intensity tails. The measured focusing intensity distribution near the design wavelength of 532 nm for the 2mm-diameter RGB-achromatic metalens is shown in Fig. S10(A). The metalens is chromatic near the design wavelength as can be seen from the focal shift. Within a 20nm bandwidth, the light is well focused. The measured focusing intensity, normalized to the one at 532 nm, is also shown in Fig. S10(B). It is possible to not only realize multi-wavelength-achromatic focusing but also to alter the dispersion near the design wavelength by using a hybrid lens system that combines a metasurface with a refractive lens. 14 Additional USAF imaging results by the RGB-achromatic metalens (NA = 0.7) Figure S11. Imaging results of the USAF resolution target group No.7 elements No.5 and 6. under synthesized color illumination. Figure S11 shows additional imaging results of the elements No.5 and 6 from the group No.7 of the USAF resolution target by the RGB-achromatic metalens. The smallest linewidth is 2.2 µm. The illumination color is synthesized by mixing blue and green light (A), blue and red light (B), and green and red light (C). The imaging results under synthesized illumination prove the RGB-achromatic imaging performance of the metalens. As long as the metalens has diffraction-limited and achromatic imaging performance for primary RGB colors, it will have the same imaging performance for other synthesized colors. It implies a platform for display imaging.  Figure S12 shows the measured focal intensity line profile of the polychromatic metalens with NA = 0.3 at six design wavelengths. The Airy function profiles are also shown in dashed lines. The discrepancy between the measurements and the Airy function profiles, which shows up at the focal intensity peak tails, is due to the large linewidth of the super continuum laser. Figure S13. Measured spectra of the super continuum laser when the output light intensity peak is tuned to center at six wavelengths.

Focusing measurement of the polychromatic metalenses
The measured linewidths of the supercontinuum laser are summarized in Fig. S13. Here, the supercontinuum laser is tuned to center the peak intensity at wavelengths of 490 nm, 520 nm, 540 nm, 570 nm, 610 nm, and 650 nm, respectively. The linewidth (FWHM of the spectrum) ranges from 4.5 nm to 7 nm.  focal intensity line profile along radial direction, which has a good agreement with Airy function profiles. The measured FWHMs are shown in the main text.

Forward design vs. Inverse design
We used our 1cm-diameter RGB-achromatic metalens design as a benchmark to compare forward design methods with our inverse design method. The results are summarized in the following Table.S3 including the objective functions and simulation results of focusing efficiencies at design wavelengths. Here, ∆ is the phase difference between a design phase and a realized phase by a meta atom at ith wavelength: i.e., ∆ = ( − , 2 ) and is in the range of -π to π. The math operator denotes the remainder after division. Ti is the normalized transmission coefficient of a meta atom at i th wavelength. For a fair comparison, we use the same meta library for both methods. In a forward design, a proper meta atom is chosen from the meta atom library to be placed at each lens position according to an objective function. Each atom is optimized independently. The objective function #1 considers only phase matching conditions and shows only 1~2% focusing efficiency at three design wavelengths. This is because the transmission of a meta atom is not uniform at the design wavelengths, and a good phase matching is realized at the cost of low and non-uniform transmission. The objective functions of forward designs #2 -#4 contain not only the phase matching condition but also the transmission. They aim to decrease the phase errors and increase the transmission of meta atoms simultaneously. The corresponding results show significant improvements compared to using objective function #1. The highest optimal focusing efficiency is 22% at wavelength of 658 nm. However, the optimization results vary with the definition of the objective functions, and the optimized focusing efficiencies are not uniform at RGB wavelengths. The above optimization results reveal the challenges when using a forward design to solve a design problem that involves multiple goals and is subject to multiple constraints. In this case, the forward design results in lower and non-uniform focusing efficiencies at achromatic wavelengths due to the lack of systematic feedback loops for optimization (in another word, it is a one-way design).
In comparison, our inverse design method achieves uniform 24% focusing efficiencies at RGB wavelengths. The meta atoms are optimized during iterations according to the feedbacks of objective function evaluation that directly reflects the performance of a desired functionality. Figure S18 (A) Our current optical setup for the laser back-illuminated micro-LCD system. (B) Our proposed optical setup for a compact display system, where a waveguide-based illumination plate is employed.

Display optical setup and strategies for improvement
In our current laser back-illuminated micro-LCD setup as shown in Fig. S18(A), we use red, green, and blue laser diodes as the RGB illumination sources. The RGB laser output beams are combined and coupled into a common fiber coupler after reflection and alignment through mirrors. The light emitted out from a fiber tip is then reflected and collimated by a reflective collimation mirror. The collimated RGB laser light goes through a light diffuser and is then incident onto the display panel. The form factor of the laser back-illuminated micro-LCD can be reduced, and a possible strategy is discussed in Fig. S18(B). In this step, we propose a waveguide-based laser back-illumination plate. An array of meta scatterers is patterned on one side, and grating couplers are patterned on the other side. The red, green, and blue laser light is first coupled into the waveguide through individual grating coupler. The coupled light is then guided through the plate and interacts with meta scatterers multiple times. The meta structures will scatter the guided wave into free space, and the light diffusion angle can also be engineered by meta scatterers. The scattered RGB laser light then illuminates onto the display after passing through a light diffuser. In this proposed setup, this laser illumination plate is more compact and has more engineering freedom to realize uniform back-light illumination.